# ---
# Fudenberg, Gao & Liang 
# "How Flexible is that Functional Form? Quantifying the Restrictiveness of Theories"
# 
# Application 3: Microfinance Takeup

# "NetNonlinearRest.R"

# Research assistant: Stephanie Nam
# date: 09/24/2023

# NOTE: This R script computes the restrcitiveness of the partially linear model

# ---

library(paramtest)
library(dplyr)
library(expm)

load(file = "NetNonlinearData.RData")
sink('NetNonlinearRest.log')

# Initialize simulated take up rates
set.seed(1)
M <- 2
p_m <- data.frame(replicate(M, runif(43))) 

# Initialize parameter ranges
theta0_range <- seq(-25, 15, by = 1)
theta1_range <- c(1e-4, 1e-3, seq(0.01, 0.5, by = 0.01))
#theta1_range <- seq(0.01, 0.3, by = 0.01)

# theta0_range <- seq(-2, -1, by = 1)
# theta1_range <- seq(0.01, 0.04, by = 0.01)

# Compute powers of adjacency matrices

A1 <- networks
A2 <- lapply(networks, function(A) A%^%2)
A3 <- lapply(networks, function(A) A%^%3)
A4 <- lapply(networks, function(A) A%^%4)
A5 <- lapply(networks, function(A) A%^%5)
A6 <- lapply(networks, function(A) A%^%6)

# Define functions

calc_H <- function (theta1, A1, A2, A3, A4, A5) {
  H <- theta1 * A1 + theta1^2 * A2 + theta1^3 * A3 + theta1^4 * A4 + theta1^5 * A5 # + theta1^6 * A6 
  return(H)
}

calc_x_j <- function (theta0, H, non, lead) {
  x_j <- c()
  for (i in 1:length(non)) {
    x_j <- append(x_j, theta0 + sum(H[lead, non[i]]))
  }
  return (x_j)
}

calc_p_ij <- function(x_ij) {
  p_j <- 1 / (1 + exp(-x_ij))
  return (p_j)
}

calc_f_i <- function(t_0, t_1) {
  H_i <- mapply(calc_H, theta1 = t_1, A1, A2, A3, A4, A5)
  x_ij <- mapply(calc_x_j, theta0 = t_0, H = H_i, non = nonleaders, lead = leaders)
  f_i <- c()
  for (i in 1:43) {
    f_i <- append(f_i, 
                  sum(mapply(calc_p_ij, x_ij = x_ij)[[i]]) / length(nonleaders[[i]])
    )
  }
  return (f_i)
}

MSE_func <- function(iter, t_0, t_1, outcome) {
  fm <- lm(outcome ~ calc_f_i(t_0, t_1))
  MSE <- deviance(fm) / 43
  return(MSE)
}


#########################################################################
############################ Restrictiveness ############################ 
#########################################################################



MSE_min_M <- data.frame()
err_nl_M <- data.frame()
err_naive_M <- data.frame()
for (m in 1:M) {
  y <- p_m[,m]
  #outcome <- p_m[,m];
  MSE_min_m <- grid_search(MSE_func, params = list(t_0 = theta0_range, t_1 = theta1_range), outcome = y)
  MSE_min_new <- cbind(MSE_min_m[[2]][which.min(unlist(MSE_min_m[[1]])),],
                    error = MSE_min_m[[1]][[which.min(unlist(MSE_min_m[[1]]))]])
  MSE_min_M <- rbind(MSE_min_M, MSE_min_new)
  
  reg_opt_m <- lm(p_m[,m] ~ calc_f_i(MSE_min_new$t_0, MSE_min_new$t_1))
  err_nl_m <- deviance(reg_opt_m)/43
  err_nl_M <- rbind(err_nl_M, err_nl_m)
  
  err_naive_m <- var(p_m[,m])*42/43
  err_naive_M <- rbind(err_naive_M, err_naive_m)
}


print("range of minimzers")
range_t0 <- c(min(MSE_min_M$t_0), max(MSE_min_M$t_0))
range_t0 
range_t1 <- c(min(MSE_min_M$t_1), max(MSE_min_M$t_1))
range_t1

mean_err_nl_rest <- colMeans(err_nl_M)
mean_err_naive_rest <- colMeans(err_naive_M)

rest_nl <- mean_err_nl_rest/mean_err_naive_rest # Restrictiveness


var_err_nl_rest <- apply(err_nl_M, 2, var) * (M-1)/M
var_err_naive_rest <- apply(err_naive_M, 2, var) * (M-1)/M
cov_err_rest<- cov(err_nl_M, err_naive_M) * (M-1)/M

avar_rest_nl <- 1/ (mean_err_naive_rest^2) * (var_err_nl_rest - 2 * rest_nl * cov_err_rest + rest_nl^2 * var_err_naive_rest)
se_rest_nl <- sqrt(avar_rest_nl /M) # Standard error for restrictiveness

print("Restrictiveness")
REST_nonlinear <- data.frame(rest_nl, se_rest_nl)
REST_nonlinear 

save.image(file = "NetNonlinearRestResults.RData")

sink()
